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I. 



INTRODUCTION 



With modern computers it is possible to simulate turbu- 
lent motion in three dimensions. As a consequence of theo- 
retical studies as well as experimental results, the 
turbulent processes at work in the three-dimensional case 
are reasonably well understood. In recent years, however, 
studies have also been conducted on two-dimensional turbu- 
lence as an analogue to the large scale motions of the atmos- 
phere. The behavior of the two-dimensional turbulence 
parameters is just beginning to be discovered. 

The principal distinguishing feature of two-dimensional 
turbulence is the tendency of its energy -containing eddies 
to grow indefinitely in scale. In any finite difference 
model of two-dimensional turbulent flow, therefore, the 
energy transfer between large scales will be computed more 
or less correctly. At the scale of the grid interval, how- 
ever, aliasing errors can transfer energy erroneously to all 
of the larger scales destroying the sense of the calculation. 
This problem can be alleviated by the removing of energy 
from scales near the grid scale. The energy thus "dissipated** 
is that which would cascade to smaller scales in real flow. 

A common procedure in numerical modeling is to para- 
meterize the sub-grid scale dissipation in the form of a 
diffusion process with an assigned coefficient of **eddy 
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diffusion**, A, as in V • (AVC) , where C is the fluid property 
being dissipated. 

The objective of the present work is to compare the 
behavior of certain properties generated by different forms 
of the eddy diffusion term in two-dimensional turbulent flow. 
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II. THE MATHEMATICAL MODEL 



The principal equation of the model is the barotropic 
vorticity equation including a space- and time-dependent forc- 
ing function and a laminar-type viscosity, written as 



8 ? 



I + + 3 = S + V • (AV^) 

u d X 



( 1 ) 



The wind velocity is associated with the stream function 
by the following equation: 

( 2 ) 



w = ii M + i,- M 
' 3y 9x 



The vorticity is the laplacian of the stream function, 

^ (3) 



In all of the numerical experiments undertaken in this 
study, the regime of integration is taken to be a rectangle 
with cyclic boundary condition in the x direction and free 
slip boundary conditions on the northern and southern bound- 
aries. (See Appendix A). 

The numerical integration of Equation (1) and (3) was 
performed by first replacing them with a pair of finite dif- 
ference equations to be solved for the variables prescribed 
on a cartesian grid. (See Appendix B) . 

The numerical space-differencing method used on the non- 
linear terms of Equation (1) is the second-order scheme 
developed by Arakawa (1966) which conserves kinetic energy 
and enstrophy within the nonlinear terms. The laplacians 
in Equations (1) and (3) were approximated by the methods 
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discussed in Appendix B. Time integration was performed by 
the explicit centered leapfrog scheme, except for the first 
step and restarts every 50 time steps by the iterative Euler- 
backward scheme. (See Appendix D) . 

The forcing function S used in these experiments was a 
crude approximation for the baroclinic driving forces in the 
atmosphere. S was negatively related to the stream function 
and was scaled to the magnitude of the rate of change of 
vorticity. (See Appendix A). The source term was set to 
zero at mid-time in each experiment conducted to study the 
decay properties of each form of dissipation. 

The initial stream field was defined to have components 
in wave numbers 2, 4, and 6 and of such magnitude as to 

simulate the 500-mb patterns. (See Appendix A). 

The different forms of the dissipation are discussed in 
the following chapter. 
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III. FORMS OF DISSIPATION 



In Equation (1), the last term on the right hand side 
represents a laminar- type viscosity or dissipation and is 
written as V • (AVC) • The experiments assumed three differ- 
ent forms of A which represents the coefficient of ”eddy 
diffusion • ” 

In the first form, A was taken to be the constant .5 x 
6 2 — 1 

10 m sec . This value was used to produce an equilibrium 
solution of a one-wave stream function field within the time 
frame of these experiments. 

The second form of A was taken to be a function of the 
wind field as formulated by Smagorinsky (1963) for three- 
dimensional turbulence: 

Ag = |d^| (Ax)^ (4) 



where is the deformation on the grid Ax and is a con- 
stant . 

The deformation is defined as: 

|I>»I - + Ds 



where 



D 



T 



8u 8 V 

8x 8y 



and 

D = ^ iu 

S 8x ^ 8y 

The third and final form of A was that proposed by Leith 
(1969) for two-dimensional turbulence, namely, 
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( 5 ) 



= Y2 Ivc1 (Ax)^ 

Here is a constant and |Vc| is the scalar magnitude 
of the gradient of vorticity. 

The finite difference forms of the above equations may 
be found in Appendix B, 
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IV, DISCUSSION OF RESULTS 



Seven experiments were conducted for comparison and 
evaluation of the different forms of dissipation as follows: 

(1) Equation (1) with the source and dissipation terms 
removed • 

(2) A constant coefficient for A in Equation (1). 

(3) Smagor insky ’ s form of A (Equation 4) with = ,5, 

(4) Smagor insky ’ s form of A with y^ = .2. 

(5) Smagor insky ’ s form of A with = .13. 

(6) Leith’s form of A (Equation 5) with y^ = .5. 

(7) Leith’s form of A with y^ = .2, 

In each experiment, the time-dependence of total kinetic 
energy, total enstrophy, and the total magnitude of the 
gradient of vorticity squared were calculated, along with 
the spectral distribution of energy and enstrophy. The 
effects of each dissipation form on these measurements and 
on typical wave disturbances were of prime interest. 

Comparing first the time— dependence of kinetic energy 
(Figures 2, 3, 4), it was noted that except for the zero dis- 

sipation, zero source case there was a tremendous decrease of 
kinetic energy initially. The order in which this dissipa- 
tion changed into regeneration was, by experiment numbers, 

7, 2, 3, 5, 4, 6. In general the Leith form halted the initial 
rapid dissipation of kinetic energy earlier. 

The energy in the model at mid-time was, by decreasing 
order of magnitude, 1 , 5 , 7 , 2 , 4 , 3 , 6 . These values were 
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obviously dependent upon the size of the dissipation term. 
However, the energy present at end-time was not in the same 
order as above, 1 , 5 , 7 , 4 , 3 , 2 , 6 . This indicated the A = con- 
stant form had a greater rate of dissipation than the other 
forms after the source term was taken out. An interesting 
point noted was that for both A^ and A^ , the smaller Y terms 
produced the greater rate of dissipation. This was due to 
the dependence of both A and on the strength of the flow 

Li U 

and with the smaller y the flow was stronger as will be 
discussed later. 

Looking next at the time-dependence of enstrophy (Figures 
5,6,7), in each of the experiments it was easily seen that 
enstrophy was far more conservative than energy. However, 
the enstrophy dissipated rapidly initially as did the kinetic 
energy. In the Smagorinsky and Leith forms of A there was a 
tendency for the enstrophy to stabilize faster for larger 
values of y. Leith's form caused a slight growth in enstrophy 
after the initial drop but once the source was cut off at 
mid-time enstrophy dissipated slightly. In Smagor insky ’ s 
form, enstrophy was more nearly conserved while the source 
was in Equation (1), but enstrophy had a tendency to be dissi- 
pated when the source was turned off. A constant A did not 
conserve enstrophy . 

The final time-dependent term to be considered was the 
total gradient of vorticity squared (Figures 8,9,10). This 
term was even more conservative than the two terms considered 
before. There was a tendency for equilibrium with respect 
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to jVcI^to be reached faster the greater the y terms. 

The fastest equilibrium solution was from experiment 6 with 
experiment 3 following, Tke experiments in which the total 
gradient of vorticity squared was most conserved were both 
Leith forms and the Smagorinsky form with Y“ • 5 , while there 
were a number of small fluctuations in 2, The deviations 
noted in the gradient of vorticity squared terms seemed to be 
more like step functions than smooth curves after the initial 
level-off of values in all experiments other than 1, 

The next measurements to be considered were the spectral 
distributions ' of kinetic energy and enstrophy (Figures 11-18) 
The scalar energy spectrum and the scalar enstrophy spectrum 
were plotted on logarithmic scales at the initial time, 1512 
time steps, and 3024 time steps. At the 1512th time step the 
source term was taken out of Equation (1) . The 3024th time 
step was the final time step. The values of kinetic energy 
and enstrophy used in the spectral analysis were taken off 
the middle grid row, K=8, Since there were 25 grid points 
east-west, it was possible to obtain 12 wavenumbers. Com- 
ponents with wavenumbers greater than 6 are subject to 
possible distortion from aliasing errors. 

In the initial spectra of kinetic energy and enstrophy, 
the energy and enstrophy were mainly contained in the even 
wavenumbers, which was due to the form of the initial stream 
field, and evenly distributed in these wavenumbers. 

The zero dissipation, zero source experiment spectra had 
no discernible pattern other than a tendency for both energy 
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and enstrophy to be spread fairly evenly in all wavenumbers. 
Experiment 2, where A was a constant, concentrated energy 
and enstrophy in the lower even wavenumbers without any dis- 
tinctive pattern. 

In experiment numbers 3, 4, and 5, there was a tendency 
for a -3 power slope to become apparent in the plots of the 
spectral variation. It was noted in the form, as y became 
smaller the spectral distribution tended to ”bend" away from 
the -3 power slope. The dropoff at the large wavenumbers end 
of the spectrum also became more pronounced as y became smal- 
ler. Obviously, the amount of energy and enstrophy present 
in corresponding wavenumbers was increased with decreasing y* 

Experiment numbers 6 and 7 had a sawtooth appearance with 
the even numbered waves generally containing more energy than 
the odd numbered waves. The envelope of the spectral distri- 
bution did have a tendency toward the -3 power slope, with the 
strongest tendency appearing at mid-time. Enstrophy tended 
to approach the -3 power slope more nearly than kinetic energy 

For a comparison of the effects of the different dissipa- 
tion forms on the wave disturbances, the stream function field 
was analyzed at the initial time and after time steps 756, 
1512, 2268, 3024 or approximately every 5 days (Figures 19-47) 
All experiments were started with the same initial field. 

The zero dissipation, zero source experiment may be referenced 
in the figures . 

After time step 756, the stream function field was begin- 
ning to become quite zonal in character. It appeared to have 
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primarily long wave components that were dominated by wave- 
number 2. The strongest flow pattern was in the experiments 
with the smallest y the A = constant case. Low areas 

were situated in the zonal band around the middle grid row, 
K=8. High areas were along the north-south boundaries. 

At time step 1512, mid-time for all experiments, the flow 
pattern was very zonal. The low areas again were along the 
middle grid rows thereby causing easterly flow in the north 
and westerly flow in the south. The experiments with the 
most zonal flow were 7, 5 and 2. Leith’s form of dissipation 
with Y= . 5 produced the least zonal flow with a wave form more 
obvious than the other forms of dissipation. 

After time step 2268, the flow is still predominately 
zonal; however, it had become weaker and was starting to 
deform in all experiments, the most notable cellular circula- 
tion being in A with y= . 5 . 

Li 

At the end-time of the experiments, the flow had con- 
siderably weakened. Leith’s form of dissipation, with y= . 2 
and the A=constant form stream function patterns were still 
very zonal in nature, Smagor insky ’ s forms deviated from the 
zonal flow more as the y term decreased. The Smagorinsky 
form with y^=.13 had cellular circulation corresponding to 
about wavenumber 2, Leith’s form with y=,5 had the weakest 
flow pattern and the least zonal. 
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V. CONCLUSIONS 



Leith's form of dissipation was the fastest acting of 
the dissipation forms. All of the time-dependent curves 
reached their minimum value faster with . After the dis- 
sipative effect had leveled off the source term had become 
dominant and the A forms had more of an opportunity to 
grow, with the smallest y coefficients growing most rapidly. 

The spectral curves showed that A had a tendency to form an 

1-1 

envelope which approximately fit the -3 power slope. The 
variation in the constant y in Leith’s form of dissipation 
had the greatest influence on the stream field. 

Smagor insky ’ s form for dissipation was not as fast act- 
ing as Leith’s form. The Smagor insky ’ s form was also not 
as sensitive to the constant coefficient as was A^ , i.e., 

all measurements made with different y’s indicated well be- 
haved curves that related almost one to one with each other. 
For the Ag form, the spectral distributions between wave- 
numbers 1 to 6 fit the -3 power slope fairly well, with 
y^^=.2 being the best fit. The stream function field for A^ 
was generally the least zonal of the forms of dissipation 
considered, with the strength of the flow inversely dependent 
upon the constant coefficients. The more zonal the flow, the 
stronger the gradients of the stream function appeared. 

The constant form of A had the slowest reaction time in 

the time-dependence measurements as long as the source term 

was still in Equation (1). However, A produced the greatest 

u 
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variations when the source term was removed in the time- 
dependence measurements. In the spectral measurements, the 
form tended to have tke -3 power slope in the even 
numbered waves only. The stream function field produced by 
A moved toward zonal flow at about the same rates as the 
other forms; however, once the zonal flow was established, 

A form kept the zonal flow more than the other forms. 

Due to the simplicity of the model used in this study 
and the use of synthetic initial data, it is impossible to 
conclude which form of dissipation is the best in meteoro- 
logical terms. However, this author believes that the 
differences in the experiments performed are great enough 
to justify further studies of this nature. 

Additional experiments that appear desirable are the 
following : 

First the expanding of the grid to investigate more 
wavenumbers . 

Secondly, the addition of a linear drag term on the 
right hand side, namely, -KC, to simulate the surface fric- 
tion forces at work in the atmosphere. 

Thirdly, and probably the most important, changing the 
source term to better simulate the baroclinic processes of 
the real atmosphere, with the source term having an effective 
wavenumber between 7 and 8. 

It is expected that these variations to this model would 
produce a more realistic spectral distribution, and in accord- 
ance with two-dimensional turbulence theory. Hopefully, 
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new formulation would produce a stream function field 
was meteorologically meaningful. 
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APPENDIX A 



COMPLETE DESCRIPTION OF THE MODEL 



15 + + + + + + + + + + + + + + + + + + + + + + + + + 

14 H 1 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 f- 

13 + + + + + + + + + + + + + + + + + + + + + + + + + 

12 + + + + + + + + + + + + + + + + + + + + + + + + + 

11 + + + + + + + + + + + + + + + + + + + + + + + + + 

10 + + + + + + + + + + + + + + + + + + + + + + + + + 

9 + + + + + + + + + + + + + + + + + + + + + + + + + 

8 + + + + + + + + + + + + + + + + + + + + + + + + + 

7 + + + + + + + + + + + + + + + + + + + + + + + + + 

6 + + + + + + + + + + + + + + + + + + + + + + + + + 

5 + + + + + + + + + + + + + + + + + + + + + + + + + 

4 + + + + + + + + + + + + + + + + + + + + + + + + + 

3 + + + + + + + + + + + + + + + + + + + + + + + + + 

2 H 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 1 1 1 1 1 1 1 1 h 

1 + + + + + + + + + + + + + + + + + + + + + + + + + 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 



Figure 1. Grid for the model. 



North-south direction has 15 grid points indexed with K. East-west 
direction has 25 grid points indexed with J. Distance between grid 
points was 2.5 latitude or 150 nautical miles (A5c=Ay=2.5°) . The grid 
is imposed on a constant P plane located at 45°N (K=8) . 



BOUNDARY CONDITIONS 

The southern boundary was placed at K=2 and the northern at K=14. 
The stream function was set equal to zero along these boundaries, i.e. 

nj,2) = = 0 

To insure free slip on the northern and southern boundaries 

iKJ,l) = - i|r(j,3) 
and i|;(J,15)= - 

was imposed on the model. 
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Periodicity was assumed in the east -west direction and achieved 

by: 

i}r(l,K) = t(26,K) 

and 1^(25, K)= 'lf(0,K) 

INITIAL CONDITIONS 

The initial stream function was defined as: 



♦ (J,K) = - C sin (2f) sin - C sin (^) sin (^) 

- C sin (- "" ■) sin (^) 

L W 



where 

X = (J-1) Zbc 
Y = (K-3) Ay 
W = 12 Ay 
L = 2W Ax 



J = 1,2,... ,25 
K = 3,4,...,16 

(width between north-south boundaries) 
(length of grid periodicity) 

C = 1.0 X 10^ 

The source terra was related to the stream function by 

S = - S f/C 
o 



where S = 1.0 x 10 
o 



-10 
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APPENDIX B 



FINITE DIFFERENCE EQUATIONS 

Equations (1), (2), and (3) are solved numerically by introducing 
finite-differences in x, y, and t. Sequential over-relaxation was 
used to solve for the field and is discussed in Appendix C. Time 
differencing is discussed in Appendix D, 

The vorticity, or the Laplacian of ijf, was approximated by the 
second-order expression using the center and four immediately conti- 
guous points, 

^ = [t(J+l,K)+ iKJ- 1,K)+ ilf(J,K+l)+ t(J,K-l)-4t(J,K)]/(Ax)^^(Ax)^ 

The numerical space-differencing method used on the non-linear 
term on the left hand side of Equation (1) is the second-order scheme 
developed by Arakawa (1966), which conserves kinetic energy and 
enstrophy within the nonlinear term, 

JCnC) = - j (CJ) 

J(C,t) = - £(J-1,K)] [>lf(J,K+l) - t|r(J,K-l)] 

-[C(J,K+D - C(J,K-1)] [t(J+l,K) - iKJ-1,K)]| /4(Ax)^ 

+ |C(J+l,K+l)[i|;(J,K+l)-t(J+l,K)]- S(J-l,K-l)[t|r(J-l,K)- 
+ C(J+l.K-l)[iKJ+l,K)-t(J,K-l)]- C(J-1,K+1 )[iKJ,K+D- nJ-l,K)] 

/4(Ax)^ 

+ |C(J+l.K)[t(J+l,K+l)-t(J+l,K-l)]- C(J-l,K)[il;(J-l,K+l)-nJ-l,K-l)] 
- C(J,K+l)['|f(J+l,K+l)-t(J-l,K+l)]+ C(J,K-1) [t(J+l,K-l)-tKJ-l,K-l)]| 
/4(Ax)^^ 

The linear P term in Equation (1) was approximated by the first order 
difference equation 

P = P { [>lf(J+l,K) - t(J-l,K)]/2Ax| 
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The above procedure was also used to find the finite difference form 
of the wind velocity components. 

For the final term in Equation (1), an imaginary staggered grid 
had to be imposed. This half grid was needed because the finite dif- 
ferencing procedure limited the values obtainable. The form of the 
viscosity term was 



V • 
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APPENDIX C 



RELAXATION TECHNIQUE 

The equation to be solved for the /ht field from Equations (1) 
and (3) is 

^ " F(x,y) (C-1) 

where 

F(x,y) = J(C,'lf) - P + S + V • (AVQ , (C-2) 

subject to the boundary conditions, /bt = 0 on the north-south 
boundaries and periodic in the east-west direction. 

In the numerical experiments, F(J,K) was computable on all interior 
grid points and the boundary conditions were known. Therefore, values 
of Sijf/3t(J,K) were determinable from the Poisson equation 

(J,K) = F^(J,K) (C-3) 

where 

Fj^(J,K) = finite difference form of F(J,K)(Zbc)^ (C-4) 

The first step in the iterative technique of sequential over-relaxa- 
tion for the solution of Equation (C-3) was to make a first guess of 
values of Sijf/St for every grid point to satisfy Equation (C-3). In 
these experiments this first guess after one step was zero everywhere; 
in subsequent steps the value from the previous time step was used for 
the guess field. 

For every estimate of the solution, the residual is defined as 
follows : 

R (J,K) = V ^ (J,K) - Fj^(J,K) 

where the superscript denotes the iteration. 
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A new estimate was made according to the formula: 






(N+l) 



(J,K) 



I 



(N) 



(J,K) + 



(1 + g) 



R (Ax) 



where OC is the overrelaxation coefficient and for these experiments 
equaled ,65. 

The iterations were continued until |r^(J,K) | ^ € for all J,K, 

-13 2 -1 

where C was taken to be 1,0 x 10 m sec at which time the solution 
was considered to be sufficiently accurate. 
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APPENDIX D 



TIME STEP SCHEMES 

In general, centered time differences (leapfrog) were used for all 
quantities; however, the first step was a Euler -backward time step. 

The latter was also introduced every 50 time steps to selectively damp 
high frequency waves and to prevent the separation of solutions at even 
and odd time steps caused by the leapfrog time scheme, 

LEAPFROG SCHEME 

+ 2At \j,K) 

where superscripts denote the time level, 

EULER-BACKWARD 

The Euler-backward is composed of two steps: 

1. (forward step) t J,K) = + At||- (J,K) 

(t+l)* 

2. (simulated backward) ij/ ( J,K) = ijf^*'\j,K)+ (JjK) 

Time steps were 10 minutes and a total of 3024 steps (about 3 weeks) 
were taken in each experiment. 
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■Figupe 8« Total gradiont of vorticity squnr'ed vs. - 

time for the zero forcing and zero dissination 
case and for A=constant form of dissipation 
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Figure 9» Total gradient of vorticity squared vs, 
time for the form of dissipation 
wi th ganrra= . . 2 , and . 1 3 
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Figure 12. Spectral distribution of kinetic energy 
and enstrophy at ir.id-tiriG and cnd-tirr.e 
for the zero dissipation and zero forcing 
case 
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Figure 13» Spectral distribution of kinetic energy 
and cns trophy at rr.id-tiinc and end-time 
for the A=constant form of dissipation 
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Figure 14. Spectral distribution of kinetic energy 
and cnstrophy at Hiid-ti-ue and end-tir.e 
for tno A<^ forn of dissipation v/ith gsrjr.a= 

.5 
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Figure 15* Spectral distribution of kinetic energy 
and cnstrophy at rrad-tine and end-tir.e 
for the A forn of dissipation with gairjaa= 

o2 ^ 
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Figure 19. Initial stream function field • 



(All stream fields contoured from -2.0x10"^ to +2.0x10"^ 
in increments of .2xlo'^ m^ sec“^) 
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Figure 20. Stream function field at 5 days for 

zero forcing and zero dissiparion case 




Figure 21. Stream function field at 5 days for 
A=constant form of dissipation 
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Figure 22 ♦ Strean function field at 10 days for 

zero forcing and zero dissipation case 




. Stream function field at 10 days for 
A=constant form of dissipation 
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Figure 24. Streaun function field at 15 days for 

zero forcing and zero dissipation case 




Figure 25. Strcari function field at 15 days for 
A=ccnstant form of dissipation 
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Figure 26. Stream function field at 20 days for 

zero forcing and zero dissipation case 




Figure 27. Stream function field at 20 da^v'S for 
A=constant form of dissi pat ion 
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Figure 28, Stream Function Field at 8 dayc For 
Aj^ Form oF dissipation with gniri'na-.^ 
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Pi/njre 29. Sti'cv.:.': Frsicti or. Field at 5 days For Aj^ 

Form oF dissi pation with gar;::;0.- , 2 
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Figure 30 . Stream function field at 10 clays for 

form of dissipation with gaimma=.5 
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Stream function field at: ]0 days 
form of dissipation with r,a;r.::: 



0 



53 



I 







Figure 32. 



Strecira function field at 15 days for 
Aj^ form of dinoination with gamma-. 5 




Figure 33c 



Strea;n f’nction field at 15 ovr.-r- for 
Aj^ form of dir.cipation with. aamma=.2 
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Figure 34. Stream function field at 20 days for 
Aj^ form of dissipation v/ith gamrna=.5 




Figure 35 o 



Strerm function field o,t 20 C:''yz for- 
form of dicsiira t Jon with ga:.ina-,2 
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Figure 360 Stream function field at 5 days for 
Ag form of dissipation v/ith gamma=,5 




Figixre 37 o Stream function field at 5 days for 

form of dissication with rarrima=,2 
s 
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Figure 38* Stream function field at 10 days for 
Ag form of dissipation with gairjua=,5 




Figure 39, Streem function field at 10 days for 
Ag form of dissipation with gamr.a= . 2 
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Figure 40, Stream function field at 15 days for 
Ag form of dissipation with gamraa=,5 




Figure 41o Stream function field at 15 days for 
Ag form of dissipation v/itli gamHia=.2 
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Figure 42. Stream function field at 20 days for 
A_ form of dissipation with gamma=,5 




Figure 43* Stream function field at 20 days for 

form of dissipation with gam:na=.2 
s 
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Figure 44» Stream function field at 5 days for 

A form of dissipation v/ith gamma=ol3 
s 




Figure 45a Stream function field at 10 days for 

form of dissipation with gamma=.13 
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Figure 46, Stream function field at 15 days for 
A form of dissipation v/ith gamma =,13 

o 




Figure 47. Stream function field at 20 days for 
A form of dissipation with gaama=,13 
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